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Abstract 

Cells store information in DNA and in stable programs of gene expres- 
sion, which thereby implement forms of long-term cellular memory. Cells 
must also possess short-term forms of information storage, implemented post- 
translationally, to transduce and interpret external signals. CaMKII, for in- 
stance, is thought to implement a one-bit (bistable) short-term memory re- 
quired for learning at post-synaptic densities. Here we show by mathematical 
analysis that multisite protein phosphorylation, which is ubiquitous in all eu- 
karyotic signalling pathways, exhibits multistability for which the maximal 
number of steady states increases with the number of sites. If there are n 
sites, the maximal information storage capacity is at least log 2 [{n + 2)/2j 
bits. Furthermore, when substrate is in excess, enzyme saturation together 
with an alternating low/high pattern in the site-specific relative catalytic effi- 
ciencies, enriches for multistability. That is, within physiologically plausible 
ranges for parameters, multistability becomes more likely than monostability. 
We discuss the experimental challenges in pursuing these predictions and in 
determining the biological role of short-term information storage. 

Key words: Multisite protein phosphorylation / multi-bit information stor- 
age / steady state analysis / multistability 
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Introduction 

Cells are information processing machines which require various forms of informa- 
tion storage, or memory, to carry out their functions. Genetic information, encoded 
in the structure of DNA, can be stably retained over an evolutionary timescale. Cel- 
lular differentiation also requires memory — whether as an adaptive response to the 
environment by a microbe or as part of the developmental process in a multicellu- 
lar organism — because information about the differentiated state must be retained 
stably after the cue which gave rise to it has disappeared. In this case, the memory 
is implemented not in the structure of a molecule but in the dynamical behaviour 
of a gene regulatory network iflTTl . Several have now been identified ll44l [P71 l24ll . 
The signalling pathways that initiate gene expression may also play a significant 
role in implementing such memories (3] 1571 . These molecular circuits all have sim- 
ilar designs consisting of interlinked positive (or double negative) feedback loops 
and mathematical analysis shows that they exhibit bistability: they have two sta- 
ble states, corresponding to two differentiation phenotypes lfT3l . Their information 
storage capacity, measured in bits, is one. To date, neither experiment nor analysis 
has shown a multi-bit capacity. 

Cells must also possess forms of post-translational short-term memory, for pro- 
cessing external signals. This is most evident in neurons. In hippocampal CA1 cells, 
transient high-frequency (tetanic) stimulation can enhance a synapse's response to 
normal stimulation. Such "long-term potentiation" (LTP) is thought to underlie 
neuronal learning and memory |[33ll25l . LTP can persist for an hour or more after 
tetanic stimulation, in a manner independent of protein synthesis (early-phase LTP), 
while repeated tetanic stimulation results in protein-synthesis dependent synaptic 
remodelling (late-phase LTP). Early-phase LTP requires a post-translational short- 
term memory. Crick and Lisman independently suggested the reaction scheme in 
Figure [IK, in which a protein kinase autophosphorylates when activated by single- 
site phosphorylation (9] [351. Lisman's mathematical analysis showed that under 
phosphatase saturation this positive feedback scheme exhibits bistability. Subse- 
quent work implicated CaMKII as a one-bit molecular memory behind early-phase 
LTP [36J. This autophosphorylating multimeric kinase is highly concentrated in 
the post-synaptic density and exemplifies the view that proteins are computational 
elements which orchestrate cellular information processing fl5][T0l. A recent model 
which presents a synthesis of current data on LTP suggests that a multi-bit capac- 
ity may be needed, although an appropriate implementation has not yet been found 

El. " 

Other signal transduction pathways, initiated by hormones, cytokines or growth 
factors, must also process complex external signals to make appropriate decisions. 
Engineering theory shows that machines with memory can undertake more com- 
plex symbol processing than machines without memory [22]. As capabilities have 



Multi-bit information storage 



3 



increased for subjecting cells to complex signals, evidence has grown for post- 
translational memory mechanisms. Point stimulation of MCF7 cells by beads coated 
in epidermal growth factor (EGF) results in rapid all-or-none activation of EGF 
receptors throughout the plasma membrane [55 J. Mathematical and experimental 
analysis of the double negative feedback loop between EGF receptor activation and 
tyrosine phosphatase activation by reactive oxygen species shows a bistable mech- 
anism underlying this ll46ll . 

In this paper we show by mathematical analysis that multisite phosphoryla- 
tion and dephosphorylation systems, which occur ubiquitously in all eukaryotic 
signalling pathways, can exhibit many stable states and that the maximal number 
of steady states increases with the number of sites. The corresponding reaction 
scheme, which requires no overt positive feedback, is shown in Figure[Tj3. If n is the 
number of sites, the maximal information storage capacity is at least log 2 (n + 2)/2 
bits, when n is even and log 2 (n + l)/2 when n is odd. If the system is initiated 
with unphosphorylated substrate then, depending on the rate constants, it can reach 
a different steady state to when the substrate is fully phosphorylated and we give in 
the Discussion an informal argument to account for this behaviour. Multistability 
predominates over monostability within physiological ranges, provided substrate is 
in excess, the kinase and phosphatase are saturated and the site-specific relative cat- 
alytic efficiencies follow an alternating low/high pattern. Furthermore, the memory 
can be switched between stable states by modulating the activity of either kinase or 
phosphatase. Our results emerge from an analytic solution for the steady state of 
the system in Figure[[jB, without the need for any rapid equilibrium or quasi-steady 
state approximations. 

Multi-bit systems can be built from one-bit systems, as in electronics. However, 
in the absence of wires and insulation, the number of components required in vivo 
would scale with the number of bits. Synthetic biologists may hence also be inter- 
ested in a molecular device with only three components which can store several bits 
of information lfTTll54l . 

Results 

Preliminary discussion of the model 

We consider a kinase E and a phosphatase F acting distributively and sequentially 
on a substrate S with n phosphorylation sites. An enzyme acts distributively if it 
makes at most one modification (addition or removal of phosphate) in each molec- 
ular encounter, so that each phospho-form competes for the enzyme. A system 
is sequential if sites are phosphorylated in a specific order and dephosphorylated 
in the reverse order. Sequentiality reduces the number of phospho-forms from 2 n 
to n + 1 and simplifies the analytical treatment developed here. If Si denotes the 
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phospho-form with i sites phosphorylated in order, then these assumptions lead to 
the chain of enzymatic reactions in Figure [Tj3. Each enzyme acts through a standard 
biochemical mechanism, as shown in Figure [IB, along with the rate constants ap- 
propriate for mass-action kinetics [8]. ATP is assumed to be kept constant by some 
external mechanism, which is not explicitly modelled, and its effect absorbed into 
the rate constants. 

These assumptions are customary in studies of multisite phosphorylation ll23l 
[321 1481 [391 [TBI 1431 but their relevance to experiment needs to be clarified. Sev- 
eral distributive enzymes have been characterised. Both Mek phosphorylation and 
MKP3 dephosphorylation of Erk, on two sites, are distributive [fT4l [6l[58l, so that 
the Mek, MKP3, Erk system is an example of a kinase, phosphatase, substrate sys- 
tem that satisfies one of the two assumptions. Sequential kinases have also been 
characterised. For instance, GSK3, in its primed phosphorylation mode, phospho- 
rylates SXXXS repeat motifs on each serine residue in a strictly C to N order |fT9"l . 
FGFR1 has also been shown to autophosphorylate in a strictly sequential manner 
|[T6l . Although these observations suggest that cognate phosphatases may act in a 
similar way, no such phosphatase is currently known. However, unlike distributiv- 
ity, which is essential for our results, sequentiality is a mathematical convenience. 
We find that non- sequential systems also exhibit multistability (not shown). We ex- 
pect this to show the same general properties as for sequential systems, although 
the maximal number of steady states may be different. 

The model has an analytic solution for the steady state 

The reaction scheme in Figure [JjB gives rise to a dynamical system of 3n + 3 or- 
dinary differential equations which describe the time evolution of n + 1 phospho- 
forms, 5*o, ■ • • , S n ; 2n enzyme-substrate complexes, ESi for < i < n and FSj 
for < j < n; and 2 free enzymes, E and F. Since the system is closed, the 
total amounts of substrate, [S tot ], and enzymes, [E tot ], [F tot ], are conserved during 
any time evolution. The system is at steady state if production and consumption of 
each species is balanced. A steady state is stable if any small perturbation causes 
a return to the state, as for a ball in a valley; it is unstable if some small perturba- 
tion causes the system to run away, as for a ball perched on top of a hill ll2~0l . The 
system is multistable if there is more than one stable steady state having the same 
total amounts of enzymes and substrate. The last proviso is important: if the system 
is initiated with different total amounts of enzymes and substrate then, because the 
amounts are conserved, it will necessarily find different steady states. This triv- 
ial possibility must always be discounted when discussing multistability in systems 
with conserved quantities. 

We showed in previous work [[TBI that this model has an analytic solution at 
steady state, without the need for rapid equilibrium or Michaelis-Menten or approx- 
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imations as customarily used [|39ll43ll49l . We briefly recall the argument. Let [— ] 
denote concentration in any steady state. Balancing production and consumption 
for enzyme-substrate complexes, we find that 



lxs,]=V§^, (i) 



where either X = E and Q<i<n or X = F and < i < n. Here, Kf denotes 
the site-specific Michaelis-Menten constant, which, using the notation in FigureQjB, 
is given by 

= £ + (2) 



Now consider the enzymatic chain in Figure [T)B. If < i < n, the net flux of 
substrate into Si from the left is always equal to the net flux out of to the right. 
For all the phospho-forms to be at steady state, it is necessary and sufficient that the 
net flux into Si from the left must equal the net flux out of Si to the right. Since 
there is never any net flux into So from the left or net flux out of S n to the right (for 
which sequentially is essential), it is necessary and sufficient that all the net fluxes 
are 0. Equivalently, each individual loop in the chain is at steady state. It follows 
that 

[gj+ii _ x a ™ 

[Si] W () 

where A; is the site-specific relative catalytic efficiency 




Applying © repeatedly, we see that 

lS l+1 ] = [S ]\ \i---\ l (^J . (5) 

It follows from ([I]) and © that if the system is at steady state then all 3n + 3 species 
concentrations are determined by [So], [E] and [F]. Conversely, if [So], [E] and 
[F] are given arbitrary positive values and the remaining species concentrations are 
defined by (OQ) and © then it can be readily shown that the system is at steady state. 
Equations (0Q) and © provide an analytic solution for any steady state of the system 
in Figure [l)B. 



Multiple steady states exist 



As explained above, multistability means the existence of two or more stable steady 
states having the same total amounts of substrate and enzymes. Equations (0Q) and 
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© enable these total amounts to be calculated in terms of [So], [E] and [F\. We 
formalise this in a function $, whose properties determine whether or not the system 
is multistable. To construct $, we need to introduce three polynomial functions of 

u = [E]/[F]: 

n 

<M U ) = A oAi • ■ ■ Ai_iu l 

M u ) = 2^ jte u ( 6 ) 

i=0 n i 

w) = 2s — jFf — u ■ 

These functions have been chosen so that, using (OQ) and ©, the total amount of 
substrate is given by (omitting the arguments of the <p functions for clarity), 

[Stot] = [So] + • • • + [S n ] + [ES ] + ■■■ + [ESn-t] + [FSi] + ■■■ + [FS n ] 

= [S ] (0! + [E!\<l> 2 + [F]<h) , 
and, in a similar way, the total amounts of enzymes are 

[E tot ] = [E] (1 + [S o ]02) 
[F tot ] = [F] (1 + [S o ]0 3 ) • 

Since [Stot] is under the control of the experimenter, while [So] is determined by the 
dynamics of the system, it is preferable to work with [S to t] instead of [So], which 
we can do by using the equation for [S tot \. We can then rewrite the equations for 
[E tot ] and [F tot ] in the form of a 2 x 2 function, $, 



$i([E],[F]) = [E]il + 
$ 2 ([£],[F]) = [F](l + 



[Stot] 



' h + [E]cf> 2 + [F]cp 3 

[Stot]fo x ' ' 



0! + [E]<f> 2 + [F]cj) 3/ 



such that $!([£], [F]) = [Et t] and $ 2 ([E], [F]) = [F tot ]. [S tot ] has now become 
part of the definition of $. The system is multistable if, and only if, $ is many-to- 
one. In other words, if there are two or more pairs ( [E] , [F] ) whose $ values are the 
same. 

Suppose that rate constants are determined and the total amounts of substrate 
and enzymes are chosen, [E tot ] = A, [F tot ] = B, [S tot ] = C. To determine whether 
or not the system is multistable, it suffices to solve the pair of equations 

$!([£], [F]) = A, $ 2 ([E],[F])=B (8) 

simultaneously for [E] and [F]. This can be done numerically as described in Ma- 
terials and Methods. The solutions give all the steady states of the system for which 
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[E tot ] = A, [F to t] = B and [S to t] = C. The curves defined by ([8]) may also be 
plotted in the ([E], [F]) plane where their intersections show the steady states. Fig- 
ure^ gives an example with four sites. The corresponding rate constants in Fig- 
ure |2K appear physiologically plausible, given our current limited understanding of 
site-specific rate constants. The curves for [E tot ] = 2.8 uM, [F tot ] = 2.8 uM and 
[Stot] = 10 fjM have five intersections, giving five steady states. A separate anal- 
ysis shows that three are stable and two unstable, as indicated. These stable states 
have widely different mixtures of the phospho-forms, as shown in Figure |2jZ!. We 
simulated the corresponding dynamical system and found that unphosphorylated 
substrate reached the steady state with low [E] and high [F], fully phosphorylated 
substrate reached the state with high [E] and low [F] (the "outer" states) and a 
suitable mixture of phospho-forms reached the inner state, as shown in Figure [2jD- 
These behaviours were characteristic of the multistable systems we simulated and 
provide a method for detecting multistability experimentally. We give an informal 
explanation for the outer states in the Discussion. 

Bistability was first shown for n = 2 in [|39l . It was later claimed that no more 
than two stable states occur when n > 2 11431 . This is incorrect, as we have just 
shown. 

A simplified solution exists when substrate is in excess 

$ gives an exact solution in two dimensions for the steady states of a 3n + 3- 
dimensional dynamical system. However, numerical solution of © is computa- 
tionally expensive. It can take up to thirty seconds to find all the steady states for 
a system with four sites, making it difficult to explore the conditions under which 
multistability arises. We found by exploration that multistability occurs when sub- 
strate is in excess so we considered what happens when either enzymes or substrate 
are in excess. If enzymes are in excess, enzyme- substrate complexes are negligible 
in comparison to [E tot ] and [F tot ], Hence, [E tot ] » [E] and [F tot ] « [F], $ is one- 
to-one and the system is monostable [fT8l . If substrate is in excess, then the total 
amounts of enzyme-substrate complexes may be considered negligible in compari- 
son to [Stot] ■ Hence, we may write, approximately, 



[Stot] = [S ] + ---+[Sn] = [S ]Mu), 



where u = [E]/[F]. We can then rewrite © to get 




(9) 
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For given [E tot ], [F to t] and [S to t}, the ([E], [F]) pairs which are solutions of © are 
the steady states of the system, to within the approximation. Dividing the first 
equation by the second, and setting [E tot ]/[F tot ] = w, we see that 

+ [S t ot}h( u )) = u(4>i( u ) + [Stot}4>2(u)) , 

and so, rearranging this, 

P{u) = (u- w)Mu) + [Stot}(u<p 2 (u) - wfa{u)) = . (10) 

Since (f>i(u), 2 {u) and (j>z{u) are all polynomial functions of u, P(u) is a polyno- 
mial function of u, whose degree is n + 1. 

For each ([E], [F]) pair which is a solution to ©, u = [E]/[F] is a positive 
solution of P(u) = 0. It can be checked that the converse is also true. Hence, so- 
lutions of the approximate system © correspond exactly to positive roots of P(u). 
Suppose that 

P(u) = a n+ iu n+1 + a n u n H h a x u + a . (11) 

The coefficients a ; may be calculated from (TTOb : 



a n+ i = A ---A n _i, a = — w and, for < i < n, 
Oi+i = Ao - ■ ■ Aj_i 



;i - kw) + is tot ) - ^ 



i+l, 



(12) 



Polynomial root finding is computationally fast and we will use this to search 
for steady states. We conducted tests and chose [Stot] /[Ftot] > 5 as the limit for 
these searches. In this range, the average normalised difference between the solu- 
tion values reported by $ and by P(u) is at most 0.23, as shown in Figure [3] and 
explained further in Materials and Methods. The frequency of potential miscount- 
ing of steady states by P{u) is 0.2% (7/3385). We considered these rates acceptable 
for the random searches below. 



The information storage capacity is at least log 2 |_(n + 2)/2j bits 

A polynomial of degree n + 1 has at most n + 1 roots [|20l . However, only positive 
roots are relevant for us. Descartes' Rule of Signs [2J states that the number of sign 
changes in the coefficients of P(u) exceeds the number of its positive roots by a 
non-negative even integer. We know from (fT2l that a n+ i > and a < 0. Hence, 
if n is odd, there can be at most n sign changes, while if n is even, there can be at 
most n + 1 sign changes: 

3 sign changes 5 sign changes 

n = 3: + - + ' n = 4: H h - H — ' • 
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Accordingly, if n is odd, the maximum number of steady states is n, while if n 
is even, the maximum is n + 1. These bounds are attained because we can show 
that any polynomial like CCD) for which a n+ i > and a < can be obtained 
by arbitrary choice of [S to t] > and appropriate choice of Kf, Kf, Aj, w = 
[Etot]/[Ftot] all positive, in (IT2|) . In particular, this can be done in such a way as to 
ensure that the approximation to the exact system ([8]) is as close as required. The 
details are given in Materials and Methods. 

Suppose then that n is odd and oti, • • • a n are any n distinct positive numbers. 
The polynomial (u — ati)(u — a 2 ) •••(« — a n )(u + 1) has degree n+1 and satisfies 
o-n+i > and a < 0. Similarly, if n is even and a a , ■ • • , a n+i are any n + 1 distinct 
positive numbers, the polynomial (it — oc\) • • • (it — a n+ i) has degree n + 1 and also 
satisfies these conditions. Hence, not only can we find rate constants for which the 
above upper bounds are attained, we can also ensure that the values of u = [E]/[F] 
at the steady states are any arbitrary pre-assigned distinct positive numbers. It is 
possible that, outside the range of approximation, the system has more steady states 
than positive roots of P(u). 11561 . following on from our results, have used singular 
perturbation theory to show that there are not more than 2n steady states. However, 
we conjecture that the bounds established in this paper always hold. 

On the basis of separate tests for stability, as discussed in Materials and Meth- 
ods, we concluded that the number of stable steady states is |_(n + 2)/2j. Since 
the information storage capacity of a system with k stable states is log 2 k bits, the 
maximal information storage capacity is at least log 2 |_(^ + 2) /2J bits. Multisite 
phosphorylation and dephosphorylation systems are capable of multi-bit informa- 
tion storage whose maximum capacity increases with the number of sites. 

An alternating low/high pattern of Aj enriches for multistability 

Under what conditions on rate constants and amounts does multistability occur and 
are these physiologically plausible? As just seen, the first question is related to 
when a polynomial has many positive roots. We found this to be mathematically 
intractable, as explained in Materials and Methods. Indeed, only probabilistic an- 
swers have been found to this general class of questions. For instance, if the co- 
efficients of (fTTI) are chosen randomly from the standard normal distribution, the 
average number of real roots (ie: without restriction on the sign) is given by the Kac 
integral formula, which is approximated by 21og(n + 1)/% [fT2l . Proportionately, 
very few of the roots of a random polynomial are real; for a random polynomial of 
degree 100, the average number of real roots is only 3.56 flU. This suggests that 
high multistability, while mathematically possible, is exceedingly rare. However, it 
still leaves open the possibility that some bias in the coefficients can enrich for it. 

According to the Rule of Signs, the number of positive roots of P(u) can only 
reach its maximum value of n + 1 when the number of sign changes in the coeffi- 
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cients is as high as possible. The sign of di+i, as given by (fT2l . is the net result of 
two additive terms, each of which may be positive or negative. We will re-interpret 
these terms in the Discussion but in the special case when Kf = K[ +1 , the sign of 
<2j is determined solely by (1 — \w). Hence, for maximum sign changes, the A, 
should satisfy an alternating low/high pattern (assuming n even): 



We found this pattern in many examples with high multistability, like the system in 
Figure [2] It is not equivalent to the alternating sign condition but has the merit of 
only involving one of the parameters. 

We find that (fT3l) enriches for multistability. We take a probabilistic approach 
to demonstrating this, in the light of the mathematical results mentioned above. For 
each even n from 2 to 12 we generated 100,000 systems as follows. We chose 
\og 10 (K* in nM) randomly from the uniform distribution on [—1, 2] and log 10 \ 
randomly from the uniform distribution on [—2, 2]. We set [Stat] = 1000 nM, forc- 
ing the enzymes into saturation, and [E to t] = [F to t] = 200 nM, ensuring that sub- 
strate was in excess. We found the distribution of steady states in Figure UK where 
monostability remains more likely than multistability up to n = 12 and five steady 
states do not appear until n = 6. We then repeated the calculation with log 10 \ 
uniform on [—2, 0] for n even and on [0, 2] for n odd, following the alternating 
low/high pattern described by (fT3l) . with w — 1. The distribution shifted to that 
in Figure @}3 in which multistability is now more likely than monostability as soon 
as n > 2, the frequency of five steady states is increased and becomes non-zero 
for n — 4. Saturation plays an important role here. We took [S tot ] = 10 uM 
and [Etot] = [Ftot] = 2 uM and found that monostability is now overwhelmingly 
more likely and that the bias in Aj has much less effect (data not shown). Hence, 
within physiologically plausible ranges, substrate excess, saturated enzymes and an 
alternating low/high pattern in the relative catalytic efficiencies enriches for multi- 
stability. 

Modulating enzyme activity leads to hysteresis 

If a multisite protein phosphorylation system acts as a memory device, it is unlikely 
to be regulated in vivo by altering its initial condition. It is more plausible that 
the activity of one of the enzymes will be modulated. We simulated the dynamical 
system in Figure [21 taking it through a cycle in [E tot ] by changing free kinase a 
small amount and letting the system relax back to a steady state after each perturba- 
tion. We found hysteresis, as shown in Figure [5K- As [E to t] is increased the system 
reaches a bifurcation point EOl where it jumps abruptly to a higher branch; when 
[Etot] is then reduced, the system remains on the higher branch beyond the bifurca- 
tion point, until jumping down to a lower branch at a lower value of [E tot ] . [E tot ] can 



Ao < — j Ai > — , • ■ ■ , A 

w w 



1 

■n-l > — 
W 



(13) 
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therefore be cycled and the system switched between the outer states in Figure [2K- 
Modulation of the enzymes can rewrite the memory. This provides another method 
for detecting multistability experimentally, which is more feasible than altering rate 
constants to show hysteresis. 

Surprisingly, systems with fewer steady states can show more complex hystere- 
sis. When [Stat] is reduced to 5 fjM the system in Figure |2jB becomes bistable with 
only three steady states. However, a similar cycle in [E tot ] produces the double hys- 
teresis in Figure 0B, showing that the system finds three stable states even though 
there is only a narrow window for [E tot ] in which three stable states exist simul- 
taneously. The potential for it, however, affects the complexity of hysteresis. We 
found a similar effect in the approach to steady state (data not shown). When a 
system is close in parameter space to regions of higher multistability, these nearby 
stable states can exert a complex influence on the dynamics. When there is merely 
the potential for higher multistability, as, for instance, when n is large, the dynamic 
and hysteretic behaviour of a system may reflect that complexity, even though the 
number of steady states in the actual system is low. 

Discussion 

Summary 

We have shown that a system with three molecular components, a kinase, a phos- 
phatase and a substrate with n phosphorylation sites, can exhibit multiple stable 
steady states and thereby function as a multi-bit post-translational cellular memory. 
The maximum information capacity increases with increasing numbers of sites and 
is at least log 2 [(n + 2)/2j bits. The conditions on rate constants for multistability 
to exist are mathematically intractable but, when substrate is in excess, enzyme sat- 
uration together with an alternating low/high pattern in the site- specific relative cat- 
alytic efficiencies enriches for multistability. That is, when rate constants are taken 
within physiological ranges and randomly sampled as specified above, multistabil- 
ity becomes more likely than monostability as soon as n > 2. The different states 
of the memory can be selected by modulating the activity of one of the enzymes. 
Even if a system has low multistability relative to the maximum, its dynamic and 
hysteretic behaviour can show the influence of nearby regions of parameter space 
with higher multistability. Our results suggest two methods for detecting multista- 
bility: different mixtures of phospho-forms can pick out different steady states — in 
particular, unphosphorylated substrate and fully phosphorylated substrate can pick 
out the outer steady states — while enzyme cycling can show hysteresis. 

While these results have been framed for protein phosphorylation and dephos- 
phorylation systems, they are potentially applicable to any reversible modification, 
such as protein ubiquitination or histone methylation [031 [301, that follows a similar 
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scheme to Figure [T^. However, much less is known about multisite effects in such 
systems. 

Multistability through kinetic trapping 

We can provide some intuition at to why the three conditions of substrate excess, 
enzyme saturation and low/high pattern of Aj give rise to two outer steady states. By 
"outer", we mean those steady states which have minimum or maximum [J5]/[.F] 
value; all other steady states are "inner". Unlike the steady-state analysis presented 
above, the argument given here follows the dynamics of the system from a given 
initial condition. In contrast to the steady state, the dynamics does not have an 
analytic solution, hence our argument is an informal one. 

Suppose that a multisite system has substrate in excess over enzymes and that 
the total amount of substrate saturates both kinase and phosphatase at each site. 
These are two of the three conditions. Let us start the system in state So with all the 
substrate unphosphorylated. Since E is saturated by So, the rate of production of 
Si will immediately reach a near maximal value, which will remain nearly constant 
as long as So continues to saturate E. As Si is produced, it will become available 
to both E, to produce S2 , and F, to produce So- However, the former reaction 
will be negligible because So, being in excess, will have sequestered free enzyme 
away from Si. The latter reaction, however, will proceed, as F is unoccupied. 
What happens next depends on the relative behaviour of E and F acting in the loop 
between So and Si. Let us assume that both enzymes work approximately according 
to the Michaelis-Menten rate law and recall (HI that these take the form 

K$ + [So] ^ W + [Si] ■ (U) 
Finally, consider a third condition: suppose that the rate curve for F lies entirely 
above that for E, as shown in Figure [6] We will interpret this in terms of the 
low/high pattern below. In this arrangement of the curves, the rate of production 
of So from Si by F can rapidly rise until it meets the nearly maximal rate of pro- 
duction of Si from So by E, at which point the So to Si loop will be in steady state. 
Although there might be a leak from Si to S2, this will be small, as long as So is 
in excess, and will be immediately balanced by back flow from S 2 to Si, since F 
is not sequestered. Hence, it seems plausible that the system will come to steady 
state with a substantial amount of S , a much smaller amount of Si and very little 
else. The phospho-form distribution becomes trapped at one end of the chain. Note 
that no other arrangement of the curves will give such trapping. If the same condi- 
tions are applied to the other end but reversed with respect to E and F, then fully 
phosphorylated substrate will become trapped predominantly as S n and the system 
will have at least two steady states. The two outer steady states in Figure [2]C! show 
exactly the distribution of phospho-forms suggested here. 
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The third condition requires that, first, the initial slope of the F curve at zero 
substrate exceeds that for the E curve and, second, that the maximal (asymptotic) 
value of the F curve also exceeds the maximal value for the E curve. From (fl4)) 
these correspond to 



c [E tot ] c([F tot ] E F 

< — tTf and c o [Ftot\ <c l [t i 



tot 



Kg K? 
respectively, which may be rewritten as 

1 - X w > and - ^ > , 

where, as previously, w = [E tot ]/[F tot \. We see from (fT2l) that this forces the coef- 
ficient ai of P(u) to be positive. For the other end of the chain, we get the opposite 
effect, with a n < 0. In other words, we recover the outer terms of the alternating 
sign condition that we found above as a necessary condition for multi stability, from 
which ([T3l emerges as a special case. We see, furthermore, that the two additive 
terms in the expression for a i+ \ in (fT2l) can be interpreted in terms of the arrange- 
ment of the Michaelis-Menten curves for E and F acting between Si and Si+i- The 
particular arrangement in Figure |6] fixes the sign of a i+ i. 

This informal argument cannot be easily extended to the inner states. If sub- 
strate is prepared in an intermediate state of phosphorylation, Si, where < i < n, 
then both E and F become sequestered and saturated immediately. Substrate will 
accumulate as S^-i and S i+ i until one or or both of E and F become accessible to 
other phospho-forms. Which of these happens will depend on other rate constants 
like af and bf, which determine the dynamics, and not just on the ones which de- 
termine the steady state like Kf and Aj (or, alternatively, Kf and cf). Hence, there 
will be many routes through which a steady state is attained, making any further in- 
formal analysis challenging. The intractability of the mathematical conditions for 
multi stability presumably reflects this complexity. Nevertheless, the inner steady 
state in Figure |2fc has an unusual distribution, with substrate concentrated predom- 
inantly in even numbered phospho-forms, suggesting that a similar type of "kinetic 
trapping" continues to determine the phospho-form distribution. 



Emergent complexity in phosphorylation, dephosphorylation systems 

Phosphorylation and dephosphorylation are ubiquitous and fundamental regulatory 
processes, which occur in all organisms. It used to be thought that prokaryotes and 
eukaryotes used fundamentally different phosphorylation chemistries but a closer 
look has revealed a more nuanced picture. Bacteria predominantly, but not exclu- 
sively, use the two-component histidine, aspartate phospho-transfer process, while 
eukaryotes predominantly, but again not exclusively, rely on serine, threonine and 
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tyrosine phosphorylation ll26l l53l l38l . A more significant distinction between the 
two kingdoms may be the extent of multisite modification. Two-component sys- 
tems typically have a single phosphorylation site on the sensor and the response- 
regulator. A recent analysis of serine, threonine, tyrosine phosphorylation in Bacil- 
lus subtilis reveals a few proteins with five to eight phosphorylation sites ll37l 
and similar maximum numbers are emerging from further bacterial studies (Boris 
Macek, personal communication, 2007). Eukaryotic proteins, however, can be far 
more heavily phosphorylated: p53, for instance, has at least 16 sites which are 
known to have regulatory function [f2TTl . 

Many suggestions have been made to account for multisite modification: signal 
integration, complex logic, attachment points to assemble signalling complexes, 
structural change through electrostatic effects, allovalency, etc [|7l |2T1 [29]|. While 
these may all be relevant, it is still puzzling why quite so many sites are needed. 
A single substrate molecule with n sites may, in principle, occupy 2 n states (over 
4000 for p53) and a population of such molecules will exhibit a distribution of these 
phospho-forms. It is not clear how such complexity can be effectively regulated 
[fT8ll49ll . Moreover, the system of kinases, phosphatases and substrate is maintained 
far from equilibrium in vivo by a steady supply of ATP. This is a recipe for complex 
emergent behaviour, as our mathematical results suggest. The in- vitro reconstitution 
of a cyanobacterial circadian oscillator 11421 . which manifests itself as an oscillation 
in multisite phosphorylation, may be an instance of such emergent complexity but 
it has otherwise proved difficult to study experimentally. 

Experimental detection of multistability 

We argued in the Introduction that signal transduction systems may require post- 
translational information storage in order to interpret complex external signals. If 
so, neither the storage mechanism nor its functional significance may be experimen- 
tally detectable in vivo without the ability to control and manipulate the signals. 
This is clear from studies of LTP in neurons: without tetanic stimulation, or some 
other complex signal to induce LTP, there would be no memory process to observe. 
T-cell activation is another context where information processing tasks have begun 
to be characterised on the basis of their response to complex signals. The T cell 
receptor is capable of being both highly discriminating among antigens and highly 
sensitive to small amounts of antigen and can accomplish both tasks quickly, a feat 
which requires an intricate mixture of kinetic proof-reading and feedback flU. It 
would not be a surprise to find short-term memory requirements in this kind of im- 
munological synapse as well. While it is technically more difficult to create and 
control signals from growth factors, cytokines or hormones, the use of microfluidic 
devices is bringing about a substantial improvement in such experimental capabili- 
ties flU. 
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The major technical obstacle in vivo, however, is the need for single-cell res- 
olution. If multistability exists, different cells in a population may be in different 
steady states and a population average could smear out the very signals that are 
being sought lfT5ll3TI . While single-cell sensors of kinase activity have been devel- 
oped [47 J, it remains challenging to determine phosphorylation state in individual 
cells. 

In-vitro studies are more feasible but, outside of extract systems 11281 . kinases 
and phosphatases have usually been studied separately (and the former more so 
than the latter). Steady states, in which kinase and phosphatase are opposed, have 
not been analysed, although there is no impediment to doing so. Care may be 
needed to ensure that the ATP is kept in sufficient excess and that ADP build-up 
does not compromise the reaction scheme in Figure [T^. Continuous-flow ATP re- 
generating systems, as used for in-vitro translation, may help ||5"T1 l27ll . The main 
difficulty lies in distinguishing and quantifying all 2 n phospho-forms of a substrate 
with n sites. Antibodies can be highly selective but we have found that, even for a 
well-studied substrate like Erk with only two sites, commercial antibodies against 
intermediate phospho-forms show too much cross -reactivity for accurate quantita- 
tion. Phospho-peptide mapping by chromatographic or electrophoretic separation 
has been successful for low numbers of sites Ifl4ll58l but mass spectrometry is now 
the proteomic method of choice and shows much promise for phospho-protein anal- 
ysis [38, 52 1 . In collaboration with Hanno Steen, we are developing methods for 
resolving and quantifying all 2 n phospho-forms using a combination of iso-electric 
focussing, HPLC and mass spectrometry. If multistability in multisite phosphory- 
lation can be detected in vitro, it seems likely that nature will have exploited it in 
vivo. 

Materials and methods 

Numerical solution of $ 

If Kf, Kf and \, L are specified and [E to t] = A [F to t] = B and [St D t] = C are chosen, then 
([8]> is solved numerically in Matlab (The Math Works, Natick, MA, USA) as follows. We 
first calculate $ on a grid in the ([E], [F]) plane and use contourc on the output to deter- 
mine the sets of points satisfying [F]) = A (the [E tot ] curve) and ^([E], [F]) = B 
(the [Ftot] curve). Contourc interpolates to find these "isolines". They provide the visual 
plots in which the steady states appear at the intersections of the curves, as in Figure 2A of 
the paper. For automated searches we use a 120 x 120 grid, where log 10 of each coordinate 
is equally spaced in [—6, 6]. For manual inspection at finer resolution we use a 1200 x 1200 
grid. We then calculate the steady states via f solve, which uses an iterative nonlinear 
search starting from a specified initial condition. We separately calculate the derivatives 
of <E> (the Jacobian) and provide that to f solve to speed up the search. An appropriate 
choice of initial conditions is essential for both speed and accuracy. We found that points 



Multi-bit information storage 



16 



lying on either the [Etot] curve or the [F to t] curve provided good initial conditions, while 
other points sometimes caused f solve to diverge or return an error. We used the [E tot ] 
curve for the set of initial conditions. We first chose three points on the [E to t] curve, one 
each at either extreme of [I?]/[_F] value and the third in the middle. If, for each of these 
initial conditions, f solve returns a solution and the solutions agree to within a specified 
tolerance (usually 10~ 4 ) in each coordinate, we return that solution as the unique steady 
state of the system. If any of these conditions fails, we take every other point lying on the 
[E t ot] curve and run f solve on all of them. We count the resulting solutions as distinct if 
they differ by more than the tolerance in any coordinate. The distinct solutions are returned 
as the steady states. This protocol was fine-tuned from numerical experiments to provide 
a reasonable balance between speed and accuracy, using the visual plot and the numerical 
calculation to cross-check each other. It can still take up to 30 seconds to find all the steady 
states for a system with four sites. 

Stability of steady states 

A dynamical system is defined by a system of ordinary differential equations, dx/dt = 
f{x), where x € R m and / : R m — > E m . The Jacobian matrix, J, is given by Jij = 
dfi I dxj . According to standard theory, the stability of a steady state is determined by the 
eigenvalues of the Jacobian evaluated at the state [20]. If all the eigenvalues have negative 
real part, the state is stable; if not, it is unstable. We computed the Jacobian symbolically 
in terms of the rate constants af ,bf , cf and the steady-state species concentrations [Y] . 
For a given steady state defined by [Stot], [E], [F], we computed all the steady-state species 
concentrations using (Q~|) and ([5]>, as described above, and substituted these values into the 
symbolic Jacobian along with the rate constants. We then calculated the eigenvalues using 
Matlab's eig function. Because the total amounts of substrate and enzymes are conserved 
we ignored the three resulting zero eigenvalues in determining the stability of a steady state. 
We found that the other eigenvalues depended on all the rate constants and not just on Kf- 
and Aj, which determine the steady state. 

In tests of stability we found that if the steady states are ordered by increasing [E] / [F] , 
unstable states typically occur between stable ones, so that typically there are (n + 2) /2 
stable states if n is even and (n + l)/2 stable states if n is odd. Both cases are covered by 
|_(n + 2)/2j, where \_x\ denotes the greatest integer not greater than x. 

Approximation of $ by P(u) 

To assess quantitatively how close P(u) = is to the exact steady state solution provided 
by <I>, we proceeded as follows with n = 4. We chose Kf and Kf randomly from the 
uniform distribution on [1, 1000] nM and log 10 A, randomly from the uniform distribution 
on [-3, 3]. We set [E tot ] = [F tot ] and chose log 10 [E tot ] and log 10 [Stat] randomly from the 
uniform distribution on [0,4], corresponding to a concentration range of [1 — 10000] nM. 
We generated 10,000 such systems, for which we solved both and P(u) for the steady 
states. We found 108 systems for which the number of steady states differed between $ and 
P(u). We first set those aside but analyse them further below. For the remaining systems, 
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we calculated for each steady state coming from $ and listed them in order of 

increasing [2£]/[.F]: s\ < S2 < ■ • • < where k is the number of steady states. (We 
found k = 1 and k = 3 only, with no k = 5.) We matched these with the ordered list 
of positive solutions of P(u) = 0, d\ < ci2 < ■ ■ ■ < ctfc. We measured the discrepancy 
between the exact solution coming from <i> and the approximate solution coming from P(u) 
by calculating the average normalised difference, 



Figure [3j\ shows that for nearly 80% of the randomly chosen systems, the approxima- 
tion is good to within a < 0.1, irrespective of the values of [Stot] an d [Etot]. Figure [3j3 
shows that the approximation gets steadily better as [Stot] /[Etot] increases from 1. We took 
[Stot] /[Etot] > 5 as our cut-off. In this range, a < 0.23. 

We then considered the 108 omitted systems for which $ and P(u) differed in the 
number of roots found. A histogram of these is plotted against log 10 [Stot] /[Etot] on the 
bottom of Figure [3}3. We found 52 miscounted systems for which [Stot]/[Etot] > 5. We 
examined each of these by hand and determined, on a conservative basis, that 45 of them 
were caused by numerical errors in <£. That is, when these systems were re-computed with 
finer tolerances and a denser set of initial conditions, the number of steady states was found 
to converge and to agree with those obtained from P(u). The remaining 7 systems were 
adjudged to be possible errors arising from using P(u) as an approximation for Since 
there were 3385 systems for which [Stot] /[Etot] > 5, this gives a miscounting rate for P{u) 
of 0.2%. 

Any polynomial can be P{u) 

Suppose given any polynomial 



with real coefficients such that A n+ \ > and Aq < 0. We claim that for appropriate choice 
of Kf, Kf , \ and w = [Et t]/[Ftot], as well as [Stat] chosen arbitrarily, all positive, the 
corresponding P(u) polynomial defined by (fl2l coincides with Q(u). We show this by 
induction. 

Note first that the term in square brackets in (fT2l can be rewritten as 



Start by choosing [Stot] > arbitrarily. Choose w = —Aq > 0. For < i < n — 2, choose 
Kf, KK_i and Aj inductively so that 



Q(u) = A n+l u n+1 + A n u n + ■■■ + Am + A 



(16) 
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as follows. (When i = 0, the induction starts with A\ = B\ but the argument below is 
identical.) If B i+ i = 0, take Aj = 1/w and choose Kf = Kf +l > arbitrarily. If 
B i+ i > 0, choose Kf > so that 



which may also always be done. If B{ + \ < then Kf, 1 may be chosen arbitrarily and the 
left hand side of (fT71 will always be positive. Hence, Kf +l and Aj can always be chosen 
positive so that (fT71) is satisfied. 

By following this inductive procedure for < i < n — 2 we have chosen [Stot], w, Kf 
for < i < n — 2, Kf for 1 < i < n — 1 and Aj for < i < n — 2 all positive. With these 
choices we have satisfied (fT2l) for all coefficients Ai such that < i < n. Now consider the 
last two coefficients A n and A n+ i. Choose A n _i = A„ + i/(Ao • • • A n _2) > 0, so that (fT2)) 
is satisfied for A n+ i. Now choose iir„_i and Kf such that 



as follows. The right hand side consists of terms like A n , which are given, or terms that 
have been previously determined. Let a = X n -iw > 0. We have to find x, y > such that 



Since a > 0, this can always be done for any c, thereby satisfying (fT2l) for A n . This 
completes the induction. 

Numerical solution of P(u) = 

We used Matlab's roots function, which is extremely fast and accurate. For n up to 12 
sites, ~ 6000 polynomials per second can be solved, giving a substantial improvement over 
numerical solution of <3?. 

Intractability of conditions for positive roots of P{u) 

Sturm's Theorem ll50l provides an algorithm for calculating the number of positive roots of 
a polynomial. We implemented this in the following Mathematica code (Wolfram Research, 
Champaign, IL, USA): 




which may always be done. Now choose Kf +l and \ so that 




(17) 




x — ay = c . 



n 




AM 

fk[u-} 



d u fo[u) 

— PolynomialRemainder[/fc_ 2 ['u], fk-i[u],u] ■ 
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Since the degree reduces by one with each remainder, the fk[u] must become constant for 
some k < n. Let v(x) be the number of sign changes in the list fo[x], • • • , f n [x]. Sturm's 
Theorem states that if a < b are not roots of /q then the number of distinct roots of /o in 
[a, b] equals v(a) — v(b). We applied this to the general polynomial a 3 u 3 + a 2 u 2 + aiu + a , 
corresponding to the case of just two sites, where we assumed that 03 > and ao < in 
accordance with (fl"2l ). We took the range to be [0, 00), using the fact that, for sufficiently 
large b, v(b) is determined by the leading coefficients of /o, • • • , fn- We found that the 
general polynomial has 3 positive roots if, and only if, the following conditions collectively 
hold: 

ai ao 2 1 do \ 

01 >0, - a o + -^<°> 9(" 3ai + ^) >0 

Qa%(a\a2 — 4afa3 + 18aoaia2a3 — ao(4af> + 27aQa^)) 
4(a| — 3aia3) 2 

These show that the region in the space of coefficients which gives rise to the maximum 
number of positive roots is highly complex. The complexity increases extremely rapidly 
with n. For n = 4 the conditions are so unwieldy that even Mathematica cannot easily 
compute them. We concluded that the question of when multistability occurs for a given set 
of rate constants is mathematically intractable. 



Model simulation 

We used the little b computational infrastructure (Mallavarapu, Thomson, Ullian & 
Gunawardena, submitted, 2007) to generate differential equation models. Little bis 
a modular programming language in which models can be specified at a biological level 
of description and compiled into Matlab code, which can then be simulated. The system 
in Figure QJ} was described inalittle b program, which was then instantiated for the 
required number of sites, making it unnecessary to write new Matlab code for different 
values of n. Little b is freely available as open source software from littleb.org and 
vcp.med.harvard.edu. For simulations we used Matlab's odel5s solver with absolute 
tolerance of 10~ 35 . 
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Figure 1: Reaction schemes for cellular memory. A Lisman scheme [35J in which 
kinase K autophosphorylates in its active state K*. B Scheme considered here, with 
no explicit positive feedback. Substrate S with n phosphorylation sites is phospho- 
rylated by kinase E and dephosphorylated by phosphatase F. Both enzymes act 
distributively and cooperate to maintain a sequential order. Si denotes the phospho- 
form with i sites phosphorylated in sequence. Phospho-forms So, • • • , S n -i have 
access to E and phospho-forms Si, - ■ ■ , S n have access to F through similar reac- 
tion schemes, with the reversible formation of enzyme- substrate complexes, ESi or 
FSj, respectively, and irreversible formation of product. With mass-action kinetics, 
each reaction has the indicated rate constant: (a for "association"; b for "break-up"; 
c for "catalysis"). ATP is assumed held constant and its effect absorbed into the rate 
constants. 
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Figure 2: Multiple steady states for a system with four phosphorylation sites. A Rate 
constants, rounded to three decimal places. Kf, Kf and A; are needed to determine 
the steady states; the other rate constants are needed to determine stability. B Plots 
of $i([£], [F]) = 2.8 fjM and Q 2 ([E], [F]) = 2.8 fM with [S tot ] = 10 fM, showing 
five steady states at the intersections. Filled squares are stable and labelled 1 (red), 
2 (black) and 3 (blue); open squares unstable. Log scales on both axes. C Bar chart 
of [So], • • • , [S4] in /jM, labelled by subscript on the horizontal axis, for each of the 
three stable states, as previously labelled. D Time courses of S 4 reaching its three 
different stable values from initial conditions [So] = a[S iot ], [S4] = (1 — a) [Stot] 
and [X] = for all other species X, with a chosen randomly in [0, 1], obtained by 
model simulation. Log scales on both axes. 
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Figure 3: Approximation of $ by P(u). 10000 random systems were generated, 
as described in the text, and solved using $ and P(u) = 0. For those which gave 
the same number of steady states, the discrepancy between the solutions was mea- 
sured using cr, as described in the text. A Histogram of log 10 o values. B The top 
shows a scatter plot of log 10 a on the left vertical axis against log 10 [S to t}/[E tot \. 
The bottom shows the number of systems which gave different numbers of steady 
states for $ and P(u), using the lower part of the right vertical axis, binned against 
log 10 [S tot }/[E 

tot • 
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Figure 4: Frequency distributions of steady states for randomly chosen systems with 
n = 2 to 12 sites, as described in the text. A Aj is chosen uniformly from site to 
site. B \i is biased to be low for even i and high for odd i. Vertical scales show 
frequency of occurrence of 1 (black), 3 (red) and 5 (blue) steady states, for 100,000 
systems for each n. 
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Figure 5: Hysteresis for the system in Figure[2j\. [E tot ] is taking in a cycle shown by 
the grey arrows, as described in the text. Note log scales on both axes. A The system 
with [E tot ] = [Ftot] = 2.8 /jM and [Stat] = 10 uM, as in Figure [2£, has 3 stable 
states. The vertical line shows its position and the numbers 1-2 mark its positions 
on the two branches of the cycle and label the corresponding steady states on the 
inserted ([E], [F]) plot from Figure |2j\- B The same system with [Stat] = 5 uM has 
only two stable states but occupies three during the cycle. The numbers 1-4 mark 
positions along the cycle — for 1, [E tot ] = 1.41 fiM; for 2 and 4, [E tot ] = 2.04 jjM; 
for 3, [E t ot] — 5-37 fjM — and also the corresponding steady states on the ([E], [F]) 
curve inserts. Changing [E tot ] alters the [E tot ] curve but keeps the [F tot ] curve fixed. 
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Figure 6: Intuitive explanation for the outer steady states. Michaelis-Menten rate 
curves are shown for E producing S\ from S (lower curve) and F producing S 
from Si (upper curve). The system is started with substrate unphosphorylated in 
state So, so that [So] is high, as shown, and Si is produced at a nearly maximal 
and constant rate, indicated by the dotted line. F is unoccupied and [Si] rapidly 
increases (grey arrow) until the rates balance, indicated by the open square, giving 
rise to a steady state. This arrangement of the curves leads to ai > in P(u), as 
explained in the text. 



